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ABSTRACT 

Magnetic fields are known to be dynamically important in the interstellar medium 
of our own Galaxy, and they are ubiquitously observed in diffuse gas in the halos 
of galaxies and galaxy clusters. Yet, magnetic fields have typically been neglected in 
studies of the formation of galaxies, leaving their global influence on galaxy formation 
largely unclear. Here we extend our magnetohydrodynamics (MHD) implementation 
in the moving-mesh code Arepo to cosmological problems which include radiative 
cooling and the formation of stars. In particular, we replace our previously employed 
divergence cleaning approach with a Powell 8-wave scheme, which turns out to be 
significantly more stable, even in very dynamic environments. We verify the improved 
accuracy through simulations of the magneto-rotational instability in accretion disks, 
which reproduce the correct linear growth rate of the instability. Using this new MHD 
code, we simulate the formation of isolated disk galaxies similar to the Milky Way 
using idealized initial conditions with and without magnetic fields. We find that the 
magnetic field strength is quickly amplified in the initial central starburst and the 
differential rotation of the forming disk, eventually reaching a saturation value. At 
this point, the magnetic field pressure in the interstellar medium becomes comparable 
to the thermal pressure, and a further efficient growth of the magnetic field strength 
is prevented. The additional pressure component leads to a lower star formation rate 
at late times compared to simulations without magnetic fields, and induces changes in 
the spiral arm structures of the gas disk. In addition, we observe highly magnetized 
fountain-like outflows from the disk. These results are robust with numerical resolution 
and are largely independent of the initial magnetic seed field strength assumed in the 
initial conditions, as the amplification process is rapid and self-regulated. Our findings 
suggest an important influence of magnetic fields on galaxy formation and evolution, 
cautioning against their neglect in theoretical models of structure formation. 
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Magnetic fields are ubiquitous in the Universe and appear 
to be present whenever ionized gas is involved. They are ob- 
served on vastly different scales, from the interiors of stars 
to interstellar gas, from haloes of galaxies to galaxy clusters 
(see, e.g., Parker 1979; Beck et al. 1996; Kulsrud 1999; Car- 
illi & Taylor 2002; Govoni & Feretti 2004). There have even 
been claims of an indirect detection of magnetic fields in cos- 
mic voids (Neronov & Vovk 2010), although the robustness 
of this result is debated (Broderick et al. 2012). 



1 INTRODUCTION 



and strength of magnetic fields in galaxies determines the 
propagation of cosmic rays (e.g. Strong & Moskalenko 1998; 
Narayan & Medvedev 2001), which are another crucial in- 
gredient influencing the gas dynamics of galaxies. In fact, the 
contribution of cosmic rays to the total pressure in the ISM 
is often assumed to be in equipartition with the magnetic 
field (Zweibel & Heiles 1997; Cox 2005), and the coupled 
dynamics of both components may be responsible for some 
of the galactic outflows (Ipavich 1975). 



In galaxies, magnetic fields are suspected to be partic- 
ularly important as here the magnetic pressure in the inter- 
stellar medium (ISM) becomes comparable to the thermal 
pressure. Magnetic fields may hence be dynamically relevant 
for the evolution of galaxies (Beck 2009), and for the reg- 
ulation of their star formation. In addition, the structure 



Magnetic field strengths have been measured for a num- 
ber of galaxies using different methods, including Zeeman 
splitting in maser emission (Robishaw et al. 2008) and radio 
polarization measurements (see, e.g., Beck 2007). These ob- 
servations provide us with quite detailed information about 
the magnetic field in the Milky Way (see, e.g., Jansson & 
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Farrar 2012a,b) and a few other nearby galaxies (e.g. Beck 
2007). 

However, it is still not well understood how galactic 
magnetic fields are originally generated and amplified, and 
which processes are most relevant for their evolution (for a 
recent review, see Kulsrud & Zweibel 2008). Weak initial 
magnetic fields could have a cosmological origin, or are cre- 
ated by Biermann batteries. A further amplification of these 
fields can then proceed through structure formation flows 
(e.g. Dolag et al. 1999), a Galactic Dynamo (e.g. Hanasz 
et al. 2004), or turbulent amplification (e.g. Arshakian et al. 
2009). Since the amplification of the magnetic field in galax- 
ies is strongly interwoven with the dynamical state of the 
gas, numerical simulations offer one of the best possibilities 
to study the complex magnetic field amplification and to 
clarify its role in regulating star formation on the scale of 
whole galaxies. 

In clusters of galaxies, magnetic fields have been in- 
cluded in a number of cosmological simulations of cluster 
growth. Dolag et al. (1999, 2002) have shown that the com- 
plex shear flows and large adiabatic compression involved 
in building up the intracluster gas lead to a sizable am- 
plification of magnetic fields, producing an end state that 
is largely insensitive to the initial seed field strength and 
configuration, and qualitatively matches observed Faraday 
rotation maps. 

So far, only few high-resolution simulations of galaxy 
formation have attempted to include magnetic fields, largely 
owing to the technical challenges involved. Wang & Abel 
(2009) studied the evolution of magnetic fields in the for- 
mation of an isolated dwarf galaxy with radiative cooling, 
but without star formation. Dubois & Teyssier (2010) stud- 
ied magnetic fields in dwarf galaxies including cooling and 
star formation, highlighting the dispersal of magnetic fields 
into the intergalactic medium (IGM) by supernova driven 
wind. These groups used a finite volume discretization of 
magnetohydrodynamics (MHD), but there have also been 
attempts to employ particle-based representations of MHD. 
In particular, Kotarba et al. (2009) compared two differ- 
ent MHD smoothed particle hydrodynamics (SPH) codes in 
simulations of magnetic field amplification in spiral galaxies. 
Kotarba et al. (2010, 2011) applied MHD-SPH techniques in 
studies of the magnetic field amplification in galaxy merg- 
ers, and Beck et al. (2012) employed SPH to simulate the 
evolution of the magnetic field in a Milky Way-like halo em- 
bedded in a cosmological environment. However, the latter 
study focussed on predictions for the diffuse gas in the galac- 
tic halo and did not discuss the properties of the magnetic 
field on the scale of the gaseous disk. 

In this paper, we aim to investigate the role of magnetic 
fields in the evolution of Milky Way-sized galaxies, includ- 
ing the regime where magnetic forces become dynamically 
important in the dense ISM gas. We shall focus on ideal- 
ized isolated galaxy models in this work, in preparation for 
future cosmological simulations of galaxy formation that in- 
clude magnetic fields. Another goal of the present study lies 
in describing and demonstrating the technical improvements 
in our updated MHD solver, which now reaches an accuracy 
that allows it to solve difficult MHD hydrodynamical prob- 
lems such as field amplification through the magnetorota- 
tional instability. 

Accounting for ideal magnetohydrodynamics in the con- 



text of galaxy formation simulations is technically chal- 
lenging, primarily because of the well-known difficulties to 
maintain the V • B = constraint in simple discretization 
schemes, an issue which is especially severe in the light of the 
very large dynamic range and spatial adaptively that is re- 
quired in galaxy simulations. One of the best approaches to 
address this problem lies in constrained transport schemes 
coupled to adaptive mesh refinement (AMR). However, this 
technique also comes with disadvantages such as compara- 
tively large advection errors for supersonic bulk flows, some- 
thing quite common in cosmic structure formation. The 
moving-mesh code Arepo (Springel 2010) is a new approach 
that overcomes this limitation, but requires a more compli- 
cated unstructured mesh geometry. Recently, we presented 
a first implementation of MHD in the Arepo code (Pakmor 
et al. 2011), which combines the accuracy of mesh-based 
techniques for hydrodynamics with the automatic adaptiv- 
ity and geometric flexibility of particle based techniques such 
as SPH. Here, we shall first begin by presenting an impor- 
tant improvement of our previous implementation of ideal 
MHD, as well as describing the required extensions for cos- 
mological integration in comoving coordinates. We will also 
show validation tests of the technique that demonstrate that 
the V • B errors are kept under control and that complicated 
flows such as strong field amplification in shear flows are cal- 
culated correctly. This method is then applied to simulations 
of isolated galaxy models, with and without magnetic fields. 

This paper is structured as follows. In Section 2, we de- 
scribe the technical implementation of our improved MHD 
solver. In Section 3 we apply our code to simulate a mag- 
netized accretion disk and compare the linear growth of the 
magneto-rotational instability to the analytical solution and 
the results of previous work in the literature. In Section 4, 
we briefly summarize the implementation of cooling, star for- 
mation, and associated feedback in our runs, and we specify 
our initial conditions. Section 5 is devoted to a discussion 
of our primary results for high-resolution simulations of iso- 
lated Milky Way-sized dark matter halos with and without 
magnetic fields. Finally, we give a summary and discussion 
of our findings in Section 6. 



2 MHD IMPLEMENTATION 

2.1 Cosmological ideal magnetohydrodynamics 

The equations of ideal MHD in physical coordinates r are 
given by 

^+V r .(pv) = 0, (!) 
^ + V r • (pvv T + p t ot - BB T ) = 0, (2) 

-^+V r -[v(£ + ptot)-B(vB)] = 0, (3) 

|+V r .(Bv T -vB T )=0. (4) 

The time derivatives are here defined at constant physical 
position r, and the spatial derivatives in V r are defined with 
respect to r. Here, p to t = p gas + §B 2 is the total gas pres- 
sure, and E — pu t h + §/>v 2 + |B 2 is the total energy per 
unit volume, with u t h denoting the thermal energy per unit 
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mass, p, v and B represent the local gas density, velocity 
and magnetic field strength, respectively. 

In addition, the magnetic field B has to fulfill the diver- 
gence constraint, VB = 0. Analytically, this constraint will 
automatically be met at all times if it is fulfilled by the ini- 
tial magnetic field. Discretization errors, however, can lead 
to a non-vanishing divergence of the magnetic field in nu- 
merical simulations. We will discuss this issue and how we 
deal with it in Section 2.3. 



2.2 Extending the MHD equations to comoving 
coordinates 

For definiteness, we briefly discuss the formulation of the 
MHD equations adopted in our code for expanding space. 
In cosmological simulations, it is convenient and common 
practice to use spatial coordinates comoving with the ex- 
pansion of the Universe, such that only gas motions rela- 
tive to the local standard of rest need to be computed. As 
usual, we parametrize the global expansion of the Universe 
with the time-dependent scale-factor a(t). Besides introduc- 
ing comoving coordinates x, we also define other 'comov- 
ing' variables (denoted with a 'c' index) for velocity, density, 
pressure, and magnetic field, as follows: 



r — ax, 

u = v — d x, 

p = p c a~ 3 , 
p = p c a~ 3 , 
B B c a 2 . 



(5) 
(6) 
(7) 

(8) 
(9) 



Here v = r is the physical velocity, and u the peculiar veloc- 
ity. With the exception of the magnetic field, these replace- 
ments are the standard definitions, but different choices have 
been used in the literature to define B c . The definition we 
adopt here has the advantage of avoiding a source term in 
the induction equation, as shown below. 

The transformed equations in comoving coordinates 
based on these replacements are as follows: 

dp c . 1, 



Q . + -V x • (p c u) = 0, 
ot a 



dt a 



p c UU +Ptot,< 



B c B r 



-Pcu, 



^ + i V x • [U (E c + Ptot,c) - B c (U B c ) 



a ( 2 o 

a \ 2a 



dt 



+ 



-V x • (b c u' 

a V 



7 uB,') 



0. 



(10) 



(11) 
(12) 



Now the time derivatives are defined at constant comoving 
position x and the spatial derivatives in V x are defined with 
respect to x. We also introduced the total comoving pressure 
Ptot,c = Pgas + and the total comoving energy density 

per unit volume, E c = p c u t h + |p c u 2 + ^B 2 . 

To get rid of most of the source terms we finally intro- 
duce scaled variables for the momentum and energy, viz. 

w — cm, (13) 
8 = a 2 E c . (14) 



Substituting these new variables into the cosmological MHD 
equations leads to the conservative system we actually solve 
in our code: 



1 

ot a 



dp, 
~dt 
dp c w 
~dT 

d8 
dt 
OB 

dt 



(p c u) - 0, 



+ V x • PcUU +p to 



+ aV x • u(E c +p t ot,c) 

iv x • (: 

a V 



1 



+ 



B c u J 



B c B r J 



B c (u B c 



0. 



o, 



2 Bc ' 



(15) 
(16) 
(17) 
(18) 



Note that in these equations only one source term remains 
as part of the energy equation. 

The above equations have a similar form as the ordinary 
MHD equations (1) in fixed coordinates, and in fact reduce 
to them for a(t) — 1, as expected. Note that we can use 
a normal Riemann solver to calculate the fluxes without 
the need for cosmological adjustments, using the following 
procedure: 

(i) We reconstruct the comoving primitive variables for a 
cell-interface in its rest-frame, using the MUSCL-Hancock 
approach as described in Springel (2010). 

(ii) We scale the magnetic field as B^ = Bc/y^- 

(iii) We then calculate the fluxes across the moving in- 
terface from the reconstructed primitive variables and B' c 
using the HLLD riemann solver (Miyoshi & Kusano 2005), 
as described in Pakmor et al. (2011). 

(iv) We finally revert the scaling of the magnetic field 
fluxes using Fb = Y' B x y/a, and multiply all other fluxes 
with the appropriate powers of the scale factor a. 

The source term in the energy equation (15) is treated in 
a Strang-split fashion by applying two half-timesteps before 
and after evolving the homogenous system by one step, sim- 
ilar as the treatment of the gravitational source terms. 



2.3 The divergence constraint 

Although the induction equation (18) analytically fulfills the 
divergence constraint of an initially divergence free magnetic 
field at all times, this property is usually lost for the dis- 
cretized version of the induction equation. However, a mag- 
netic field which develops a non-zero divergence can lead to 
severe artifacts in numerical simulations. In particular, di- 
vergence errors can induce a huge local amplification of the 
magnetic and velocity fields, rendering both unphysical. 

An elegant solution to this problem is the constraint 
transport approach (Evans & Hawley 1988). This scheme 
discretizes magnetic and electric fields such that the se- 
quence of their update steps manifestly guarantees a van- 
ishing divergence of the magnetic field at cell centres, to 
machine precision. Unfortunately, however, this method is 
at present limited to Cartesian grids, and it is conceptually 
unclear whether it can be adapted to unstructured moving 
meshes at all. 

Instead of trying to completely avoid divergence errors 
by construction, as in constraint transport, we therefore at- 
tempt to keep them as small as possible, so that they do 
not affect our results. In Pakmor et al. (2011) we imple- 
mented the Dedner divergence cleaning approach (Dedner 
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Figure 1. Orszag-Tank vortex at t = 0.5 for the Dedner divergence cleaning (top row) and the Powell scheme (bottom row). Shown 
from left to right are density, pressure, magnetic energy, and the absolute value of the relative divergence error as defined in Eqn. (23). 
Both simulations were run with 600 x 600 cells. 



et al. 2002) in the moving-mesh code AREPO and showed 
that it is usually able to keep the divergence error at a small 
level. In the Dedner scheme, a local divergence error is both 
advected away from the place where it originates and also 
damped at the same time. However, as shown in Pakmor 
et al. (2011), a relatively restrictive timestep criterion is re- 
quired to use the Dedner scheme together with individual 
timesteps for the cells, and even then it is difficult to guar- 
antee stability in very dynamic environments, encumbering 
applications to cosmic structure formation. 

To remedy these issues, we now adopt the Powell ap- 
proach for divergence control (Powell et al. 1999). In this 
scheme, additional source terms are introduced into the mo- 
mentum equation, induction equation and energy equation: 



+ - V x • (pell) = 0, 

ot a 

<9p c W / T 

h Vx • I p c uu +p t ot,< 



dS 

dt 



+ aV x 



(V x -B c ) B c , 



U (E c +p t ot,c 



1 



B c (U B c ) 



|B C 2 - i (V x • B c ) (u • B c 

-7T- + -V x B c u - uB c 
ot a V 



(19) 



(20) 



(21) 



(Vx-B c )u. (22) 



These source terms encode a passive advection of V • B/p 
with the flow and counteract further growth of local V • B 
errors. By experience, we find that this scheme is very stable 
in practical applications. A welcome advantage is also that 
it is completely local and does not require any additional 



constraints on the timesteps, even when all cells are evolved 
on individual timesteps. 

We apply the divergence source terms for all active cells 
at the same time when we calculate and apply the local 
fluxes. In each timestep, we first define the magnetic field 
perpendicular to a cell interface as the average of the per- 
pendicular magnetic fields left and right of the interface, 
B x = | (B x ,l + B x ,r). This value of B x is then used in the 
Riemann solver. We estimate the divergence of the magnetic 
field of an active Voronoi cell by 



V-B, 



(23) 



using the average normal component B x at the interfaces 
calculated as above, and add the fluxes and divergence 
source terms at the same time to our conservative variables. 
To estimate the relative size of the divergence error, we typ- 
ically also use the values for the divergence of the magnetic 
field calculated as part of this treatment of the source terms. 



2.4 Testing the divergence control 

As one important test of our new implementation of diver- 
gence control based on the Powell approach, we simulate 
the Orszag-Tang vortex problem with the same configura- 
tion as shown in Pakmor et al. (2011), using both the pre- 
viously implemented Dedner divergence cleaning and the 
Powell approach. The two simulations are shown at time 
t = 0.5 in Fig. 1, for comparison. There is excellent agree- 
ment between both results in density, pressure and magnetic 
energy, although the employed approach to limit divergence 
errors is significantly different. Note that this agreement is 
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Figure 2. Maps of the divergence error for a Orszag-Tank vortex 
at t = 0.5 for the Dedner divergence cleaning (top row) and the 
Powell scheme (bottom row). The left column shows the absolute 
value of the relative magnetic divergence error on a logarithmic 
scale, whereas the right column shows its signed value on a linear 
scale. Both simulations were run with 600 x 600 cells. 



a non-trivial outcome, given that small and seemingly in- 
nocent changes in the implementation details of each of the 
divergence control methods are typically readily punished 
by clearly visible artifacts in these physical quantities. This 
is because both approaches are intricate numerical schemes 
that just barely keep the V • B related errors under control. 

There are, however, some differences in the relative size 
of the resulting divergence errors. As is clearly visible in 
Fig. 1, the absolute value of the relative divergence error is 
larger in some regions for the Powell approach compared to 
the Dedner divergence cleaning. This is visible even better 
in Fig. 2, where maps of the divergence error are shown. 
At time t = 0.5, we find a volume- weighted average of the 
magnitude of the relative divergence error equal to 9 x 10 -4 
for the Dedner divergence cleaning, and 2.4 x 10 -2 for the 
Powell approach, respectively. If the sign of the divergence 
error is not ignored in the averaging, the mean divergence 
error is much smaller, only — 3 x 10 ~ 5 for the Dedner scheme 
and —2 x 10 -4 for the Powell method, because positive and 
negative fluctuations in V • B tend to largely cancel. 

It is important to note that for both measurement meth- 
ods the average relative divergence error is much smaller 
than unity and spatially highly localized. The spatial dis- 
tribution of the errors is very similar for both divergence 
cleaning schemes. Interestingly, in places where the absolute 
value of the relative divergence error becomes large, its sign 
alternates between neighboring cells. This suggests that the 
divergence error occurs as a locally confined problem, most 
likely originating in approximate solutions of the Riemann 
problems and in the local discretization and truncation er- 
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Figure 3. Radial magnetic field at t = 30 for a test of the 
magneto-rotational instability, after about five inner orbits. 



rors. The local errors are oscillatory in nature, and the par- 
tial cancellation among neighboring cells helps to keep any 
large-scale impact very small. Although the average diver- 
gence error is about an order of magnitude larger for the 
Powell scheme compared to the Dedner cleaning, there are 
no noticeable differences in any fluid quantities, as shown in 
Fig. 1, supporting this interpretation. This good alignment 
of the two schemes in the Orszag-Tang vortex problem is 
only lost at much later time, when the initial vortex has de- 
cayed into turbulence, and the turbulent state as a function 
of time becomes slightly different between the schemes. In 
Sec. 5.7, we will return to a discussion of the divergence error 
in our applications of the code to realistic galaxy simulation. 



3 MAGNETO-ROTATIONAL INSTABILITY 

The magneto-rotational instability (MRI) (Balbus & Haw- 
ley 1991) is one of the most important MHD phenomena in 
astrophysics, making it a critical test-problem for numeri- 
cal MHD codes. In particular, in accretion disks the MRI 
is suspected to provide the primary means for an efficient 
transport of angular momentum in the radial direction, al- 
lowing the gas to be accreted by the central compact object. 
The MRI may also facilitate the growth of magnetic fields in 
galaxies (Kitchatinov & Riidiger 2004). Interestingly, for ac- 
cretion disks it is possible to calculate the linear growth rate 
of the MRI analytically (Balbus & Hawley 1991), allowing 
for a direct quantitative check of whether a numerical MHD 
implementation can correctly account for this process. 



3.1 Test setup 

We employ the same setup as used in Flock et al. (2010), but 
simulate the full disk to avoid the need for azimuthal bound- 
ary conditions. In this setup, a thin, non self-gravitating, 
Keplerian disk rotates around a central body of unit mass, 
M = 1, implying a centripetal acceleration of g(R) — 1/R 2 
in cylindrical coordinates. The disk extends from R = 1 to 
R = 4 in radial direction and from z = —0.5 to z = 0.5 
in vertical direction. Initially, it has uniform density p — 1 
and pressure P — c^p/j, with c s = 0.1 and 7 = 5/3, and 
rotates with an azimuthal velocity of = y/l/R. The ra- 
dial and vertical velocity is given by axisymmetric random 
perturbations of amplitude ±5 x 10 -4 . 

A uniform vertical magnetic field with a strength of 
0.05513/n and n — 4 pervades the disk between the radii 
R — 2 and R — 3. We use periodic boundary conditions in 
the vertical direction and a layer of boundary cells at R = 1 
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Figure 4. Maximal amplitude of the radial magnetic field at 21 
equally spaced radii from R = 2 to R = 3 versus local orbits. The 
red and green line show the amplitude at radii R = 2 and R = 3, 
respectively. 

and R — 4. These boundary cells keep the undisturbed parts 
of the disk at their initial radial position. 

The mesh-generating points are initially arranged on 
128 circular rings in the plane of rotation and duplicated 
64 times uniformly in the vertical direction. The number of 
cells per ring is set such that cells have a mass of ~ 10 -5 
per cell, with a uniform spacing within each ring. 

For this setup, the radial and azimuthal magnetic field 
should be amplified exponentially in time according to 
B = B exp(7MRi£) (Balbus & Hawley 1991). The fastest- 
growing mode in vertical direction should be the n = 4 
mode, the dominant radial mode depends on the initial per- 
turbations. The critical (fastest possible) mode grows with 
7mri = 0.7511, with an angular frequency Q = R~ 15 as a 
function of radius. 



3.2 Results 

Figure 3 shows the radial magnetic field after about 5 inner 
orbits, corresponding to 1.8 orbits at R = 2. The MRI is 
clearly active and a radial magnetic field has emerged. The 
dominant vertical mode of the instability is close to n = 4, 
as predicted for this setup. The magnetic field is largest at 
R — 2 and becomes smaller with increasing radius (for the 
pattern, compare to Flock et al. 2010). We note that the 
"checkerboard" -instability that occurs for some implemen- 
tations tested in Flock et al. (2010) does not show up in our 
simulation. 

The growth of the amplitude of the radial magnetic field 
at different radii is shown in Fig. 4. Between about 0.5 and 
1.5 local orbits, the amplitude increases exponentially, es- 
sentially at all radii. Only the very inner and outer radii at 
R = 2 and R = 3 are slightly different, most likely due to 
boundary effects. Starting after about 1.5 local orbits, the 
growth becomes steeper at a progressively larger range of 
radii, since the magnetic field is already saturated in the 
inner parts of the disk, which leads to non-linear growth in 
regions close to the saturated part. 

We fit the growth rate at different radii selected be- 
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r 

Figure 5. Exponential growth rate of the magnetic field ampli- 
tude measured at different radii (blue line). The black line shows 
the maximum growth rate expected analytically. 

tween 1.0 and 1.5 local orbits and show the result in Fig. 5. 
The growth rate in units of the inverse orbital period ranges 
from 0.53 to 0.72 at different radii with an average growth 
rate of 0.66, in very good agreement with the results by 
Flock et al. (2010) for the same setup, which were based on 
completely different MHD implementations that included 
constraint transport. Also, there are no regions in our re- 
sults where the local growth rate exceeds the fastest possible 
growth rate of 0.75. Reassuringly, we thus conclude that our 
code is able to correctly reproduce the linear phase of the 
magneto-rotational instability. 



4 SIMULATION SETUP AND ADDITIONAL 
PHYSICS 

Our disk galaxy simulations start from slowly rotating 
spherical gas clouds embedded in equilibrium into a colli- 
sionless dark matter halo (using the setup of Jubelgas et al. 
2008). The dark matter halo is modelled with a static back- 
ground potential corresponding to a NFW density profile 
(Navarro et al. 1997) with a mass of M200 = 10 12 Mq and a 
concentration parameter of 7.2 roughly resembling the halo 
of the Milky Way. We employ a spin parameter of A = 0.05 
and a baryonic mass fraction of 0.17. In the simulations with 
magnetic fields, we typically introduce a homogeneous seed 
magnetic field of 10 ~ 9 G parallel to the x-axis. The same or 
similar seed field strengths have also been used by Wang & 
Abel (2009), Kotarba et al. (2009) and Dubois k Teyssier 
(2010). Indeed, intergalactic fields of this size are expected 
in models of IGM magnetic seeding through outflows from 
dwarf galaxies at high redshift (Kronberg et al. 1999). Nev- 
ertheless, we have also varied the strength of the seed field 
in a subset of our tests. 

Gravitational forces are calculated with a standard tree- 
method (based on Barnes & Hut 1986; Springel 2005) to ac- 
count for the self- gravity of the gas. The forces due to dark 
matter are determined from the static background poten- 
tial, neglecting a possible the adiabatic contraction of the 
halo when baryons cool and settle in the centre. We use the 
Arepo code in its pseudo-Langrangian configuration for our 
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Figure 6. Time evolution of different quantities of a 10 12 Mq galaxy with (black lines) and without magnetic field (black lines) for a 
mass resolution of 2.1 x 1O 5 M0. The top left panel shows the total baryonic mass in stars and gas (straight line) and the gas mass 
(dashed line) within a radius of 15kpc. The top right panel shows the thermal pressure (dashed line), the magnetic pressure (dotted 
line) and the total pressure (straight line). The pressure is calculated as volume- weighted average in a sphere of radius 15kpc around the 
center of the galaxy. The bottom left panel contains the volume-weighted average root mean square of the absolute value of the magnetic 
field within a radius of 15kpc. Finally, the bottom right panel shows the total star formation rate in the whole simulation. 



simulations. In this mode, the mesh-generating points, which 
define the Voronoi cells, are moved with the local fluid ve- 
locity, subject to small corrections to keep the shape of cells 
reasonably regular. In addition, we use refinement and dere- 
finement operations where needed to ensure that the mass 
of the cells always stays within a certain narrow range: if a 
cell contains more than twice or less than half of the desired 
average mass per cell, we split it into two cells or merge it 
with its neighbours, respectively. The refinement and dere- 
finement operations are carried out as described in Springel 
(2010) and Vogelsberger et al. (2012). 



The gas is allowed to cool radiatively, which eventually 
causes the rotating gas sphere to develop a strong cooling 
flow and to form a rotationally supported disk inside-out 
at the centre. For simplicity, we only include atomic cool- 
ing by helium and hydrogen, ignoring molecules or metals. 
Therefore, the gas can only cool down to a temperature of 
10 4 K. We include star formation and supernova feedback by 
means of a simple sub-resolution model (Springel & Hern- 
quist 2003). This model assumes that star formation and 
associated supernova feedback lead to a self-regulated multi- 
phase interstellar medium in which cold molecular clouds are 
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Figure 7. Resolution dependence of different quantities in simulations of a galaxy with a seed magnetic field of 10 9 G. Shown are 
simulations with a mass resolution of 2.1 x 10 5 M (black lines), 1.6 x 10 6 M (blue lines), and 1.4 x 10 7 M (red lines). The top left 
panel shows the total baryonic mass in stars and gas (straight line) and the gas mass (dashed line) within a radius of 15kpc. The top 
right panel shows the thermal pressure (dashed line), the magnetic pressure (dotted line) and the total pressure (straight line). The 
pressure is calculated as volume- weighted average in a sphere of radius 15kpc around the center of the galaxy. The bottom left panel 
contains the volume-weighted average root mean square of the absolute value of the magnetic field within a radius of 15kpc. Finally, the 
bottom right panel shows the total star formation rate in the whole simulation box. 



embedded in a tenuous hot phase roughly at pressure equi- 
librium. The pressurization of the medium due to supernova 
feedback can can be described by an effective equation of 
state. Stars are formed in a probabilistic approach out of 
the medium, with a gas consumption time-scale consistent 
with the observed Kennicutt (1998) relation for star- forming 
disk galaxies in the local Universe. Once a star is formed 
from a cell, we remove 90% of the mass, momentum and 
internal energy of the cell and create a collisionless star par- 
ticle at the position of the cell that inherits the removed 



mass and momentum. Since the cell only retains 10% of its 
mass, often a derefinement of the cell is triggered in the 
next timestep, meaning that the cell will be dissolved and 
its contents merged with its neighbors. 

We do not include any local supernova feedback in the 
form of point-like energy depositions. Such local supernova 
feedback would likely increase the amplification of the mag- 
netic field, because of the implied additional shearing mo- 
tions and the higher turbulence in the interstellar medium. 

Also, we do not change the magnetic field of a cell when 
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Figure 8. Gas density and magnetic field at t = 1 Gyr. Top two rows show projections, bottom rows slices through the center of the 
galaxy. The left column shows the density for the simulation with a mass resolution of 2.1 x 10 5 M without a magnetic field. Center 
and left column show density and magnetic energy (for the projection) and magnetic field (on the slice) for the simulation with the same 
resolution but an initial magnetic seed field of 10 — 9 G. 
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Figure 9. Gas density and magnetic field at t = 2 Gyr. The two rows on top show projections, the bottom rows slices through the center 
of the galaxy. The left column depicts the density for the simulation with a mass resolution of 2.1 x 10 5 Mq without a magnetic field. 
Center and left column show density and magnetic energy (for the projection) and magnetic field (on the slice) for the simulation with 
the same resolution but an initial magnetic field of 10 — 9 G. 



© 0000 RAS, MNRAS 000, 000-000 



|v| [kms 



Simulations of magnetic fields in isolated disk galaxies 11 

|B|LuG] 



□ 



a 



10 1 



10 2 



15 



10 



z 



-5 t 
10 

-15 



1 / 1 / y 1 > '1 

, / / / / / 

/ / f 

/ / / / 

, / / / / 

' '/ / y y ' ' ^ " 
' ' 

/ ' ' ' 7 / 
- / / J ! - / 

\ . J / t t M \\ ^-O , 

v v . ^ ^ - N / 
^ v \ ' ^ / \ / - \ \ / - 


■* - "x r > 

^ x^ v /V / - /' t 

* v v ^ </ ' V, / L 

* -xx \ \ 

y \ t \ // \ x — ^ n \ 

- X X i - — / — -J - - 


- { ^ \ x -x- / f \ - 1 *- \ 

'\ N i\ / / - 
^ / x X v - - \- • ^1 

- / - - ------ / X ^Xv 

- i - . 'J Mir- 


/ \ //I \ v K - v |\ 
- " "* M " 



-15 



-10 



-5 
x [kpc] 



10 



15 




1 r i i 1 


)" 2 


lO" 1 


10° 


10 1 10 2 




^ U -^^^ \ / ^^^_x-^^ / z\ ^^^^^ 



10 



-15 



; ^ \ ^ ^ V/ ^-^^v x x ^ 



|^^t \ X^' ~ f — ^ N ^ NN nX V 




^ . \ \ xx x ^ 



10 



15 



10 



15 -10 -5 5 


10 15 


x [kpc] 




p [lO lo M kpc- 3 ] 




■^i i i i 1 

1 i i i i 1 


r 6 10- 5 io- 4 10- 3 


io- 2 io- 



-10 



-15 



-15 -10 -5 5 
x [kpc] 



10 15 



Figure 10. Properties of outflows above and below the disk at t = 2 Gyr for the simulation with a mass resolution of 2.1 x 10 5 Mq with 
a magnetic field. Shown are the magnitude of the velocity (top left panel), the magnitude of the magnetic field (top right panel), the 
ratio of magnetic pressure to thermal pressure (bottom left panel) and the density (bottom right panel) in slices perpendicular to the 
disk. 



a star is formed, similar to Dubois & Teyssier (2010) in 
their AMR-code. This has the major advantage that the 
local structure of the magnetic field is not affected by a 
star- forming event, and in particular, the divergence of the 
magnetic field does not increase when a star is formed. Note, 
however, that this is a conservative and non-trivial assump- 
tion. For the mass resolution we can afford, each "star parti- 
cle" is a macro-particle that represents a whole stellar popu- 
lation formed in a molecular cloud. In such a cloud, magnetic 
fields can be strongly amplified by turbulence. At the same 



time, some part of the magnetic field is locked up in stars, 
whereas another part is pushed out of the cloud by ambipo- 
lar diffusion. Therefore, the field dynamics in star-forming 
clouds is complicated and in reality it may well be possi- 
ble that the magnetic field surrounding a molecular cloud is 
amplified during star formation. Similarly, but arguably less 
likely, it is possible that a molecular cloud leaves a smaller 
magnetic field behind after being depleted by star formation. 
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Figure 12. Time evolution of the volume-averaged total pres- 
sure (straight lines), thermal pressure (dashed lines) and magnetic 
pressure (dotted lines). Black, green, blue and cyan lines repre- 
sent galaxies with halo masses of 1O 12 M , 10 n M o , 1O 1O M , 
and 10 9 Mq, respectively. 

5 DISK GALAXY RESULTS 

Without radiative cooling, the initial gas spheres in our sim- 
ulations remain in hydrostatic equilibrium. However, as we 
allow the gas to cool, the pressure support in the centre is 
removed and the gas sphere collapses. Since gas elements 
carry some angular momentum, they settle into a dense, ro- 
tationally supported disk that can locally fragment and form 
stars. The result is the formation of a disk galaxy in a classic 
inside-out fashion (Fall & Efstathiou 1980). Similar numer- 
ical experiments have frequently been used to study models 
of star formation and feedback (e.g. Springel & Hernquist 
2003; Jubelgas et al. 2008). Here, we are primarily inter- 
ested in whether magnetic fields can impact the dynamics 
of this disk formation scenario. 



5.1 Global evolution 

The global evolution of two galaxies with and without a 
magnetic seed field is compared in Figure 6. For the first 
1.5Gyrs the evolution is quite similar. After only about 
200 Myrs the inner part of the gas sphere has collapsed and 
formed a small disk in which a large starburst has ensued. 
In this starburst, the star formation rate peaks just below 
100M©yr _1 , but drops to about 10 - 20M©yr _1 within a 
few hundred million years, at which point a more or less con- 
stant, only slowly declining star formation rate is reached. 

The magnetic field is strongly amplified in the initial 
collapse of the gas cloud and the first starburst, from an 
initial field strength of 10 ~ 9 G to an average field strength 
of the order of fiG in the inner parts, in agreement with ob- 
servational constraints (Beck 2007). Much of the subsequent 
amplification proceeds though an Q-type dynamo, made effi- 
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cient through the continued radial gas inflow combined with 
strong differential shear. The magnetic field saturates at the 
level of few fiG, and no significant further amplification takes 
place in our simulation in places where the magnetic field 
has saturated. Once it has been amplified to this saturation 
level, we find that it also contributes significantly to the to- 
tal pressure. Owing to the additional magnetic pressure, it 
becomes slightly harder for the gas to form stars. There- 
fore, as shown in Fig. 6, the star formation rate is lowered 
by about 30% at late times in the simulation with mag- 
netic fields, compared to the run without magnetic fields. 
At late times, after more than 1.5Gyrs, the average mag- 
netic pressure even begins to dominate the average thermal 
gas pressure in parts of the star-forming gas disk. 



5.2 Resolution dependence 

Before we proceed with further analysis, we investigate the 
convergence of our results with numerical resolution. To this 
end, we have run the simulation with magnetic field at three 
different resolutions, using mass resolutions of 2.1 x 10 5 M , 
1.6 x 10 6 M , and 1.4 x 10 7 M , respectively. The evolution 
of the galaxy for these different resolutions is compared in 
Figure 7. 

The two higher resolution runs agree fairly well with 
each other. In particular, the magnetic field saturates at 
about the same value and it takes about the same time to 
reach saturation. In the lowest resolution simulation, how- 
ever, the magnetic field is amplified somewhat faster, and 
saturates at an average value that is about a factor of three 
larger compared to the higher resolution runs. The magnetic 
pressure in this simulation also takes over the thermal pres- 
sure earlier (at about t = 1 Gyr) and the ratio of average 
magnetic pressure over average thermal pressure is signifi- 
cantly larger than unity. This finally leads to a considerably 
stronger suppression of star formation, and a late-time star 
formation rate of about one tenth of the late-time star for- 
mation rates in the high-resolution runs. 

These differences suggest that in the low resolution 
run the amplification of the magnetic field is artificially in- 
creased due to poor resolution, caused in particular by rela- 
tively large divergence errors (see Fig. 13). We note however 
that the impact of these errors on the amplification is very 
much smaller than in the SPH calculations of Kotarba et al. 
(2009), where order of magnitude differences were found for 
different approaches to control the divergence errors. In con- 
trast to our results, Dubois & Teyssier (2010) have found 
slightly smaller amplification for poor resolution, perhaps 
as a result of their accurate suppression of divergence er- 
rors. For isolated galaxies, the weak systematics we see in 
our technique for low resolution is not a serious problem, be- 
cause such models can nowadays be simulated with at least 
the medium resolution we use here. In cosmological simula- 
tions of structure formation where galaxies form hierarchi- 
cally, the reduced accuracy at low resolution is potentially a 
more serious concern. Here, the simulations always contain 
small progenitor galaxies which are initially resolved only 
with typically very few cells. Accurately capturing the dy- 
namics of these systems when magnetic fields are involved 
may be quite challenging with our techniques. 



5.3 Structure of the gas disk and its magnetic 
field 

Projections of the gas density in the disk and maps of its 
magnetic field are shown in Figures 8 and 9, after one and 
two Gyrs, respectively. After one Gyr, the inner part of the 
gas halo has collapsed to a disk galaxy and converted a large 
fraction of its gas into stars. The structure of the galaxy 
at this time is very similar for the simulations with and 
without magnetic field. Even the position and morphology 
of individual spiral arms are very similar, and there is no 
difference in the thickness of the disks. This is consistent 
with the result that the magnetic pressure is on average still 
subdominant at this time, and therefore magnetic fields are 
not expected to have a substantial effect on the dynamics of 
the gas yet. 

In the calculation with magnetic fields, the field has 
been amplified to maximum values of up to 100 /iG in the 
center of the disk and to a few /iG in its outer parts. The 
strength of the magnetic field is correlated with the spi- 
ral arms, with a significantly larger magnetic field in spiral 
arms than in between arms (consistent with Kotarba et al. 
2009). Compared to the smooth spirals arms seen in the den- 
sity field, the magnetic field strength shows a more patchy 
appearance with alternating regions of larger and smaller 
magnetic field strength. Interestingly, although we do not in- 
clude winds from supernovae, the simulation with magnetic 
fields locally shows highly magnetized outflows from the disk 
which also transport magnetic fields outwards. While overall 
the magnetic field in the disk is aligned with the spiral struc- 
ture of the disk, it also shows a lot of field reversals, similar 
as in Dubois Sz Teyssier (2010). The field strength is about 
a factor of 10 larger in the central part of the disk compared 
to its outskirts. In the central part, where the timescale on 
which differential rotation in the disk amplifies the magnetic 
field is smaller, the magnetic field is also more regular. We 
note that the structure of the simulated magnetic field is in 
good agreement with observational data of spiral galaxies 
(Beck 2007; Jansson & Farrar 2012b). 

After an elapsed time of two Gyrs, the differences be- 
tween the galaxies simulated with and without magnetic 
fields have become larger. By this time, the total gas mass in 
the disk has dropped by about a factor of two, as a result of 
star formation and the reduced supply of gas cooling out of 
the halo. The galaxy without magnetic field still shows dis- 
tinct spiral arms and also retains about the same thickness 
it had one Gyr earlier. In contrast, in the simulation with a 
magnetic field, the individual spiral arms have mostly disap- 
peared to form a smoother disk with weaker residual spiral 
patterns. The magnetic field in the disk is now very regular 
and well aligned with the residual spiral structures of the 
disk. This change can be explained in terms of the longer 
differential rotation time in the outer parts of the disk. Now 
after 2 Gyrs, differential rotation had enough time to also 
amplify the magnetic field in the outer parts to saturation 
values. Only in the very center of the disk a significant ra- 
dial component remains. Field reversals in the disk still exist 
but have become much less frequent compared to the struc- 
ture of the magnetic field at time t = 1 Gyr. In addition, 
the gas disk is now noticeably thicker compared to the disk 
without magnetic fields, as a result of the magnetic pressure 
contribution. 
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Figure 13. Time evolution of the relative divergence error of the magnetic field for the simulations with a mass resolution of 2.1 x 10 5 Mq 
(black lines), 1.6 x 10 6 Mq (blue lines), and 1.4 x 10 7 Mq (red lines). Left and right panels show the volume-weighted average divergence 
error and the volume-weighted average absolute divergence error. 



The magnetic field is primarily amplified by two pro- 
cesses in our simulations. In the first few lOOMyrs (during 
the initial starburst), adiabatic compression of the magnetic 
field dominates the amplification. As material cools, it be- 
comes denser until its density reaches the star formation 
threshold. Then it forms stars and leaves a highly amplified 
magnetic field behind. Therefore, the magnetic field strength 
is initially correlated with the local star formation and is 
highest at the center of the disk, where the compression and 
the star formation rate is highest. Later, after the total star 
formation in the disk has dropped significantly, the amplifi- 
cation of the magnetic field is dominated by shear motions in 
the differentially rotating disk combined with radial inflows. 
This leads to a pronounced inside out growth of the mag- 
netic field in the disk, because the orbital timescales increase 
for larger radii. The relentless shearing of the magnetic field 
in the disk also make it ever more regular with time. 

There are other processes which are likely important for 
the amplification of magnetic fields in galaxies but which are 
not included in our present simulations. Because we evolve 
an isolated galaxy, it does not experience any large-scale 
shearing motions caused by infalling material or mergers 
with other galaxies that could lead to an additional ampli- 
fication of the magnetic field. In addition, our subresolution 
model for the interstellar medium does not directly resolve 
interstellar turbulence, because we include supernova feed- 
back smoothly through an effective equation of state (see 
Sec. 4). Therefore, we miss a possible turbulent amplifica- 
tion of the magnetic field on small scales. 

In addition, we do not explicitly distinguish between 
neutral gas, which should not feel the magnetic field directly, 
and ionized gas, which interacts with the magnetic field. 
Instead, we neglect resistivity effects and treat all gas in 
the limit of ideal magnetohydrodynamics which will tend to 



overestimate the dynamical effect of the magnetic pressure 
on the gas in the disk. 



5.4 Outflows 

The magnetized galaxy shows significant outflows of gas up 
to a few kpc above the disk and below the disk, as shown in 
Fig. 10. These outflows are strongly magnetized and their 
velocity vector is highly correlated with the magnetic field 
lines. We note that unlike in the simulations of Dubois & 
Teyssier (2010) these outflows are not driven by supernovae. 
While supernova feedback is implicitly included in our sub- 
resolution model for the ISM, it does not create galactic 
winds; instead its effects are absorbed in the applied effective 
equation of state. 

The outflows are driven by low density, highly magne- 
tized bubbles which rise above the disk. In our simulation 
they reach typical velocities of a few hundred km s _1 , almost 
reaching escape speed. The further they rise above the disk, 
the smaller their velocity and magnetic field becomes until 
they finally fall back onto the disk again. The outflows con- 
tain about 10 8 Mq of gas at t = 2 Gyr. Qualitatively, these 
outflows are similar to galactic fountain models in which 
outflows are driven by supernova-powered bubbles (see, e.g., 
Shapiro & Field 1976; Oosterloo et al. 2007; Marinacci et al. 
2011). The quantitative properties of the outflows, however, 
have to be taken with care, since the resolution of our sim- 
ulation is already rather low in these regions. 



5.5 Varying the magnetic seed field 

To understand the effect of the magnetic seed field on the 
evolution of the model galaxies, we repeat our simulation 
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Figure 14. Relative divergence error of the magnetic field at t = lOGyr (top row) and at t = 20 Gyr (bottom row). The left column 
shows the absolute value of the relative divergence error, the right column shows the signed values. 



with the intermediate resolution (which had a mass resolu- 
tion of 1.6 x 10 6 M and a seed field of 10~ 9 G) with a larger 
and a smaller seed field of 10 -8 G and 10 -10 G, respectively. 
The evolution of pressure, magnetic field and star formation 
rate compared to the original simulation is shown in Fig. 11. 

Although the initial seed field differs by a factor of 
10 compared with the original simulation, the evolution of 
the calculation with a ten times smaller seed field is hardly 
changed. In this case, the growth of the magnetic field is de- 
layed by about lOOMyrs during the initial 300Myrs. After 
this time, when the initial starburst subsides, the magnetic 
field in both simulations converges to the same value and 



evolves basically identically, indicating that the magnetic 
field saturated to a state independent of the initial seed field. 
Because even directly after the starburst the magnetic pres- 
sure is still too small compared with the thermal pressure 
to be of dynamical relevance, the different choice of the seed 
field in this case does not affect the evolution of the galaxy 
or the evolution of the magnetic field once it is saturated. 

The behavior of the simulation with a larger magnetic 
seed field of 10 -8 G, however, is slightly different. After the 
initial starburst, the magnetic field is temporarily about 
twice as large as in the original simulation (where a seed 
field of 10 -9 G was used). From this time onwards, it stays 
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Figure 15. Relative divergence error of the magnetic field at t = 10 Gyr. The top left panel shows the error sampled on a 256 x 256 
grid. Top right, bottom left, and bottom right panel show the error from the top left panel averaged over 2 x 2, 4 x 4, and 8x8 cells. 
This averaging rapidly reduces the residual local errors, indicating that the errors are locally compensated and do not induce incorrect 
large-scale magnetic dynamics. 



roughly constant for two Gyrs at an elevated level until it 
finally decreases to the same strength as in the other two 
simulations that started with one or two order of magni- 
tude weaker seed fields. Again, since the magnetic pressure 
is much smaller than the thermal pressure during the initial 
starburst and for most of the following time, this difference 
has only a small effect on the total pressure or the star for- 
mation rate. Only at late times, just before the magnetic 
field decreases to the value found in the other simulations, 
the star formation rate is suppressed slightly more than for 



the other two calculations with smaller seed fields. We hence 
find that the saturation strength of the magnetic field is 
almost completely insensitive to the initial seed strength, 
apart from residual differences stemming from the differ- 
ent times required for amplifying the field from the initial 
strength to saturation. 
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5.6 Halos of different masses 

Another important parameter in our simulations is the mass 
of the simulated galaxy. To complement our set of simula- 
tions, we have run three more simulations with 10 5 gas cells, 
with the total mass of the halo set to 10 11 M , 10 10 Mq, and 
1O 9 M , respectively, and an initial magnetic field strength 
of 10 -9 G. The virial radii of the halos were scaled in a self- 
similar fashion. 

In Figure 12, we show the evolution of the average 
pressure in the galaxies forming in these simulations and 
compare it to our fiducial simulation with a halo mass of 
1O 12 M0. The most important difference lies in the relative 
strength of the magnetic pressure developing for galaxies of 
different mass. For the two simulations with halo masses 
of 10 12 M and 1O 11 M , the magnetic pressures reaches 
equality with the thermal pressure after about two Gyrs. 
However, for the two lower mass galaxies, the magnetic field 
never becomes dynamically important. This difference from 
self-similar behavior is to be expected because of the non- 
selfsimilarity of cooling and star formation physics. Both the 
cooling efficiency and the ability of a halo's gravitational po- 
tential well to compress the gas of the interstellar medium 
against the pressure delivered by supernovae vary strongly 
with halo mass. In the lower mass halos, the potentials be- 
come too shallow for compressing the gas to very high densi- 
ties, therefore the amplification of the magnetic field by adi- 
abatic compression of the seed field is weaker. Furthermore, 
the much more anemic and thicker gas disks developing in 
these halos reduce the amplification of the magnetic field by 
differential shear. We thus expect that magnetic fields are 
particularly important only for those halos that are also effi- 
cient sites of star formation. According to recent abundance 
matching constraints for the ACDM cosmology (Guo et al. 
2010; Moster et al. 2010), these are primarily halos in the 
mass range - io 11 - 5 " 12 - 5 M . 

In the simulations with the 10 12 M and 10 11 M galax- 
ies, outflows are present as discussed in Section 5.4, but they 
are absent in the two low mass galaxies with halo masses of 
10 9 M and 10 10 M . This directly confirms that these out- 
flows only occur once the magnetic field becomes dynami- 
cally important. We stress that these outflows are not caused 
by supernova feedback, rather they can form when the mag- 
netic field pressure dominates over the thermal pressure of 
the gas in at least parts of the disk. 



5.7 Divergence error 

As discussed in Sec. 2.3, our discretization scheme of the 
MHD equations does not manifestly maintain a vanishing di- 
vergence of the magnetic field to machine precision; rather 
a finite divergence may appear, with the scheme prevent- 
ing its further growth. Although such small divergence er- 
rors do not cause any obvious artifacts even in complex en- 
vironments (including technically dicey steps such as de- 
refinement and star formation), it is important to examine 
the properties of the divergence error in more detail. We 
note that at the very least it can alter the local jump condi- 
tions for the Riemann solver, which in principle can some- 
times lead to locally incorrect fluxes, causing larger numer- 
ical noise or possibly more severe artifacts. 

In Fig. 13, we show the evolution of the average relative 



divergence error as a function of time for different numer- 
ical resolutions in our default halo. With increasing reso- 
lution the divergence error decreases, which indicates that 
it is dominated by local noise in the magnetic field. This is 
also consistent with the average divergence error being much 
smaller (by more than a factor of 10) than the average ab- 
solute divergence error. 

The spatial distribution of the divergence error is shown 
in Figure 14. Here it can be seen directly that the sign of the 
divergence error is alternating locally on small scales. There 
are basically no coherent patches of multiple cells with only 
positive or negative divergence error. The magnitude of the 
divergence error follows the pattern of the spiral arms, which 
can be easily understood as a consequence of the divergence 
errors correlating with large gradients of the magnetic field. 
In other words, large (unresolved) gradients in the local mag- 
netic field cause large divergence errors at the same place. 
With higher resolution, the gradients are resolved better and 
the divergence errors become smaller. 

Another confirmation of the local nature of the diver- 
gence error can be inferred from Figure 15, where local spa- 
tial averages of the divergence error are shown. They quickly 
reduce the error to the sub-percent level for progressively 
larger smoothing regions. This suggests that non-zero er- 
rors in the field divergence do not introduce changes in the 
large-scale dynamics of the magnetic field. Note also that 
the sizes of the divergence errors seen in our simulations are 
smaller or at most of the same order than in cosmological 
MHD simulations carried out with the latest generation of 
SPH-MHD methods (Dolag & Stasyszyn 2009; Beck et al. 
2012; Stasyszyn et al. 2012), which have similar properties 
regarding magnetic field divergence errors. 



6 CONCLUSIONS 

In this study, we have examined the magnetic field amplifi- 
cation in simple models of disk galaxy formation. Our goal 
has been to provide a first exploration of the potential im- 
pact of magnetic fields for the regulation of star formation in 
Milky Way-sized galaxies. We have chosen isolated galaxies 
in order to allow a study of the magnetic field evolution in 
a well-controlled setting, while still confronting our numeri- 
cal scheme with all the technical challenges (such as dealing 
with star formation events and on-the-fly derefinements) it 
has to cope with in future cosmological applications. 

Another important aspect of our study lies in testing 
our new MHD implementation in the moving-mesh code 
Arepo, which we extended to cosmological simulations of 
galaxy formation that account for radiative cooling and 
star formation. In contrast to our previous method (Pak- 
mor et al. 2011), we here use an 8-wave formulation (Powell 
et al. 1999) for dealing with the divergence constraint of the 
magnetic field. This new implementation turns out to be 
much more robust, whereas our previous MHD code could 
sometimes become unstable in very dynamic environments 
unless very strict timestep contraints were imposed. 

We have shown in this paper that our new scheme pro- 
duces competitive results for test problems and is successful 
in keeping the divergence error small. We have demonstrated 
that our new code is able to simulate the magneto-rotational 
instability in disks without problems, and that the correct 
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linear growth rate is reproduced. Also, we have shown that 
the new implementation prevents any noticeable artifacts in 
simulations of disk formation. 

For our simulations of isolated galaxy formation, we ob- 
serve good convergence for the strength of the magnetic field, 
the magnetic pressure and the star formation rates, unless 
the resolution is very low. In particular, we obtain converged 
results for our two high resolution runs with mass resolu- 
tions of 2.1 x 10 5 M and 1.6 x 10 6 M , corresponding to an 
initial resolution of 8 x 10 5 gas cells and 1 x 10 5 gas cells 
in the halo, respectively. The magnetic field saturates at a 
strength of about 10 — 100 ji G in the very center of the disk 
in a Milky Way-sized halo of mass 1O 12 M0, and at a few 
/iG in the main body of the disk. We find that the star 
formation rate and structure of the gas disk is very similar 
for our runs with and without magnetic fields until about 
1 Gyr. Later, the magnetic pressure becomes comparable to 
the thermal pressure, providing additional support for the 
gas and reducing the star formation rate by about 30%. It 
also causes changes the structure of the disk, reducing the 
prominence of individual spiral arms as the additional mag- 
netic pressure produces slightly more homogeneous disks. 
Interestingly, the magnetic field causes weak outflows from 
the disks which rise several kpc above the disk before their 
material falls back. 

Simulations of halos with a reduced mass of 10 11 M 
yield qualitatively very similar results. Here the magnetic 
field is also amplified to a saturation level, allowing it to 
affect the galaxy dynamics at late times. In contrast, small 
halos with masses of 10 10 M and 10 9 M did not amplify 
the magnetic field up to equipartition with the thermal pres- 
sure, due to insufficient star formation and gas compression. 
This can be understood as a consequence of the shallow grav- 
itational potential well and the modified cooling efficiency 
in these systems. 

Our results underline the potentially very important 
role magnetic fields play in galaxy formation, especially in 
systems that have a high efficiency of star formation, as here 
the gas is compressed particularly strongly and the amplifi- 
cation through gas shearing motions and inflows tends to be 
strong, too. In such galaxies, magnetic fields of typical IGM 
strength are quickly amplified to saturation values, such that 
the initial seed field strength is of very limited importance 
and tends to be quickly forgotten. It will be very interest- 
ing to examine the influence of magnetic fields on forming 
galaxies in cosmological hydrodynamic simulations of struc- 
ture formation, which we plan to study in future work. The 
present work has demonstrated that such calculations are 
technically feasible and sufficiently accurate with our present 
MHD implementation in Arepo. 
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